---------- matrix operations ------------
-- see the XFun SVD for a more efficient invert()

multiply(A,B)[i,j] = sum(A[i,k]*B[k,j],k,
         1,count(B)) dim[count(A),count(B[1])]

transpose(A)[i,j] = A[j,i] dim
                        [count(A[1]),count(A)]

invert(A) = adjoint(A)/det(A)

adjoint(A) = transpose(cofactor(A))

cofactor(A)[i,j] = (-1)^(i+j)*
                   det(submatrix(A,i,j))

submatrix(A,k,l)[i,j] = A[i,j] when i<k and j<l,
            A[i+1,j] when ik and j<l,
            A[i+1,j+1] when ik and jl,
            A[i,  j+1] when i<k and jl dim
                        [count(A)-1,count(A)-1]
